6 休息中の凝集性に影響する要因

以下の分析では、どのような要因が群れ内の凝集性に影響しているかを検討します。

6.1 データの加工

まず、休息中の周辺個体数などについて算出して元データに結合します。

  • RG_female: 休息集団内のメス数
  • RG_female_plus1: 自身を含めた休息集団内のメス数
  • RG_est: 休息集団内の発情メス数
  • RG_**: それぞれのオスが休息集団内にいたか
  • x3m_female: 3m以内のメス数
  • x5m_female: 5m以内のメス数
  • **_3m: それぞれのオスが3m以内にいるか
## 各調査期間ごとの6歳以上のメスのID  
adult18 <- unique(female_18m %>% filter(age >= 6) %>% .$femaleID)
adult19 <- unique(female_19m %>% filter(age >= 6) %>% .$femaleID)
adult20 <- unique(female_20m %>% filter(age >= 6) %>% .$femaleID)
adult21 <- unique(female_21m %>% filter(age >= 6) %>% .$femaleID)

## 10m以内近接個個体名を記したデータを作成  
focal_prox <- focal_combined_final %>% 
  separate(x0_1m, into = str_c("x0_1m",1:11), sep = ",") %>% 
  separate(x1_3m, into = str_c("x1_3m",1:15), sep = ",") %>% 
  separate(x3_5m, into = str_c("x3_5m",1:16), sep = ",") %>% 
  separate(x5_10m, into = str_c("x5_10m", 1:9), sep = ",") %>% 
  pivot_longer(cols = x0_1m1:x5_10m9,
               names_to = "proximity",
               values_to = "ID") %>% 
  mutate(proximity = ifelse(str_detect(proximity,"x0_1m"),"x0_1m",
                            ifelse(str_detect(proximity,"x1_3m"),"x1_3m",
                                   ifelse(str_detect(proximity,"x3_5m"),"x3_5m",
                                          ifelse(str_detect(proximity,"x5_10m"),"x5_10m","NA"))))) %>% 
  filter(!is.na(ID)) 

## 休息集団に含まれるメスの数  
RG_female <- focal_prox %>% 
  filter(RG == "1") %>% 
  filter((study_period == "m18" & ID %in% adult18)|(study_period == "nm19" & ID %in% adult18)|
           (study_period == "m19" & ID %in% adult19)|
           (study_period == "m20" & ID %in% adult20)|(study_period == "nm21" & ID %in% adult20)|
           (study_period == "m21" & ID %in% adult21)|(study_period == "nm22" & ID %in% adult21)) %>% 
  left_join(female_all %>% select(date, femaleID, rs2), 
            by = c("date","ID" = "femaleID")) %>% 
  group_by(date, no_focal, time) %>% 
  summarise(RG_female = n(),
            RG_est = sum(rs2, na.rm = TRUE)) %>% 
  ungroup()

## 休息集団内に各オスがいるか否か  
RG_TY <- focal_prox %>% 
  filter(RG == "1") %>% 
  select(date, no_focal, time, ID) %>% 
  filter(ID == "TY") %>% 
  mutate(RG_TY = 1) %>% 
  select(-ID)

RG_IT <- focal_prox %>% 
  filter(RG == "1") %>% 
  select(date, no_focal, time, ID) %>% 
  filter(ID == "IT") %>% 
  mutate(RG_IT = 1) %>% 
  select(-ID)

RG_LK <- focal_prox %>% 
  filter(RG == "1") %>% 
  select(date, no_focal, time, ID) %>% 
  filter(ID == "LK") %>% 
  mutate(RG_LK = 1) %>% 
  select(-ID)

RG_KR <- focal_prox %>% 
  filter(RG == "1") %>% 
  select(date, no_focal, time, ID) %>% 
  filter(ID == "KR") %>% 
  mutate(RG_KR = 1) %>% 
  select(-ID)

## 3m、5m近接(相手個体の活動に依らず)    
focal_prox_all <- focal_prox %>% 
  mutate(ID = str_replace(ID, "\\(F\\)","")) %>% 
  mutate(ID = str_replace(ID, "\\(M\\)","")) %>% 
  mutate(ID = str_replace(ID, "\\(O\\)","")) %>% 
  mutate(ID = str_replace(ID, "\\(N\\)","")) %>% 
  mutate(ID = str_replace(ID, "\\(n\\)","")) %>% 
  mutate(ID = str_replace(ID, "\\(LT\\)","")) %>% 
  filter((study_period == "m18" & ID %in% adult18)|(study_period == "nm19" & ID %in% adult18)|
           (study_period == "m19" & ID %in% adult19)|
           (study_period == "m20" & ID %in% adult20)|(study_period == "nm21" & ID %in% adult20)|
           (study_period == "m21" & ID %in% adult21)|(study_period == "nm22" & ID %in% adult21)) 
  
x3m_female <- focal_prox_all %>% 
  filter(proximity %in% c("x0_1m","x1_3m")) %>% 
  group_by(date, no_focal, time) %>% 
  summarise(x3m_female = n()) %>% 
  ungroup() 

x5m_female <- focal_prox_all %>% 
  filter(proximity %in% c("x0_1m","x1_3m","x3_5m")) %>% 
  group_by(date, no_focal, time) %>% 
  summarise(x5m_female = n()) %>% 
  ungroup() 

## 元データに結合  
focal_combined_prox <- focal_combined_final %>% 
  left_join(RG_female, by = c("date","no_focal","time")) %>% 
  left_join(x3m_female, by = c("date","no_focal","time")) %>% 
  left_join(x5m_female, by = c("date","no_focal","time")) %>% 
  left_join(RG_TY, by = c("date","no_focal","time")) %>% 
  left_join(RG_IT, by = c("date","no_focal","time")) %>% 
  left_join(RG_LK, by = c("date","no_focal","time")) %>% 
  left_join(RG_KR, by = c("date","no_focal","time")) %>% 
  replace_na(list(RG_female = 0, RG_est = 0, x3m_female = 0, x5m_female =0, RG_TY = 0,
                  RG_IT = 0, RG_LK = 0, RG_KR = 0)) %>% 
  replace_na(list(x0_1m = "NA",x1_3m = "NA", x3_5m = "NA", x5_10m = "NA")) %>% 
  ## 自身を含んだ休息集団サイズ  
  mutate(RG_female_plus1 = RG_female + 1) %>% 
  ## オスとの近接情報を追加  
  mutate(TY_3m = ifelse(str_detect(x0_1m,"TY")|str_detect(x1_3m,"TY"),1,0),
         IT_3m = ifelse(str_detect(x0_1m,"IT")|str_detect(x1_3m,"IT"),1,0),
         KR_3m = ifelse(str_detect(x0_1m,"KR")|str_detect(x1_3m,"KR"),1,0), 
         LK_3m = ifelse(str_detect(x0_1m,"LK")|str_detect(x1_3m,"LK"),1,0)) 

続いて、その日群れ内にいた血縁個体数を計算し、結合します。

kin <- read_csv("../Data/data/others/kin.csv")

female_presence <- group_all %>% 
  select(groupID, study_period, date, Kur:Yun) %>% 
  select(-c(TY,IT, LK, KR, KM, TG)) %>% 
  pivot_longer(cols = Kur:Yun,
               names_to = "femaleID",
               values_to = "presence") %>% 
  left_join(att) %>% 
  ## 6歳以上の個体のみを抽出
  filter(age >= 6,
         presence == 1) 

focal_numkin <- focal_combined_final %>% 
  distinct(study_period, date, groupID, no_focal, subject) %>% 
  left_join(female_presence %>% select(date, groupID, femaleID),
            by = c("date", "groupID")) %>% 
  filter(subject != femaleID) %>% 
  left_join(kin, by = c("subject" = "femaleID", "femaleID" = "femaleID2")) %>% 
  mutate(kin = ifelse(kin >= 0.0625,1,0)) %>% 
  group_by(no_focal, date, subject) %>%
  summarise(no_kin = sum(kin)) %>% 
  ungroup() 

focal_combined_prox %>% 
  left_join(focal_numkin, by = c("no_focal","date","subject")) -> focal_combined_prox_b

続いて、TY、ITとの親密度(CSI)の情報も結合します。

focal_combined_prox_b %>% 
  left_join(CSI_TY_combined %>% select(subject, CSI_TY)) %>% 
  left_join(CSI_IT_combined %>% select(subject, CSI_IT)) -> focal_combined_prox_c

最後に、その日観察したオトナメスへの攻撃頻度も追加します。

### 6歳以上への攻撃のみ  
aggression_all %>% 
  mutate(femaleID = str_replace_all(femaleID, "\\?","")) %>% 
  left_join(att, by = c("femaleID", "study_period")) %>% 
  ## 確実に被攻撃個体が6歳以上であるものを抽出
  filter(age >= 6 | str_detect(femaleID,"multi|many|Multi|KunTrt|Ako,Kil")) %>% 
  group_by(date) %>% 
  summarise(no_agg = length(unique(no))) %>% 
  ungroup() %>% 
  filter(date >= "2018-10-08") -> no_agg_6yo

focal_combined_prox_c %>% 
  left_join(group_all_b %>% 
              select(groupID, date, duration)) %>% 
  left_join(no_agg_6yo, by = "date") %>% 
  mutate(rate_agg_6yo = no_agg*60/duration) %>% 
  replace_na(list(rate_agg_6yo = 0)) %>% 
  left_join(female_all %>% select(date, femaleID, rs2),
            by = c("subject" = "femaleID",
                   "date")) -> focal_combined_prox_final

作成したデータは以下の通り。

focal_combined_prox_final

6.2 休息集団サイズの分析

ここでは、TYやITの有無を含む様々な要因が個体追跡中の休息集団サイズ(正確には、休息集団内のメス数)に影響するかを検討します。

6.2.1 交尾期(非発情メス)

6.2.1.1 データの加工

まず、非発情メスについて分析を行うためのデータを作成します。連続変数は標準化します。

  • 休息ポイント数が10以上の個体追跡セッションのみを使用
  • ハドリングのポイントは除去
  • 地上休息のみを含む
cohesion_data_an <- focal_combined_prox_final %>% 
  mutate(rate_female = no_female/max_female) %>%
  filter(rs2 == 0) %>% 
  replace_na(list(hud = 0)) %>% 
  filter(hud != "1") %>% 
  filter(RG == 1 & T_G == "G") %>% 
  group_by(date, no_focal, subject, study_period, rate_agg_6yo, no_ntm, 
           no_est, no_male, no_female, TY_presence,
           no_kin, IT_presence, CSI_TY) %>% 
  summarise(sum_RG = sum(RG_female_plus1),
            n = n()) %>% 
  ungroup() %>% 
  mutate(TY_presence = as.factor(TY_presence),
         IT_presence = as.factor(IT_presence)) %>%
  filter(study_period != "m18") %>% 
  filter(n >= 10) %>% 
  mutate(ntm_std = standardize(no_ntm),
         male_std = standardize(no_male),
         est_std = standardize(no_est),
         kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         rate_agg_std = standardize(rate_agg_6yo),
         CSI_TY_std = standardize(CSI_TY),
         mean_RG = sum_RG/n)

cohesion_data_an


共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。

cohesion_data_an %>% 
  pivot_longer(cols = c(no_ntm, no_est, no_female, rate_agg_6yo, kin_std, CSI_TY),
               names_to = "type",
               values_to = "value") %>% 
  ggplot(aes(x = value, y = mean_RG)) +
  geom_point(shape = 1) +
  facet_rep_wrap(~ type,
                 scales = "free") +
  theme_bw(base_size = 12) +
  theme(aspect.ratio = 1) +
  labs(x = "", y = "Mean resting group size")


6.2.1.2 モデリング

それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。

  • 応答変数: log(平均休息集団サイズ(メスのみ))
  • 分布: Studentのt分布
  • 説明変数: TYの在不在(TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
  • ランダム切片: subject
## 2019~2021年
m_cohesion_an <- brm(log(mean_RG) ~ TY_presence + IT_presence + ntm_std + est_std + 
                       female_std + kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd"),
                                           prior(student_t(4,0,10), class = "nu"),
                                           prior(student_t(4,0,10), class = "sigma")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = cohesion_data_an,
                                 file = "model/m_cohesion_an.rds")

## 2019~2020年
cohesion_data_an_bf21 <- cohesion_data_an %>% 
  filter(study_period != "m21") %>% 
  mutate(ntm_std = standardize(no_ntm),
         male_std = standardize(no_male),
         est_std = standardize(no_est),
         kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         rate_agg_std = standardize(rate_agg_6yo),
         CSI_TY_std = standardize(CSI_TY),
         mean_RG = sum_RG/n)

m_cohesion_an_bf21 <- brm(log(mean_RG) ~ TY_presence + IT_presence + ntm_std + est_std + 
                       female_std + kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = cohesion_data_an_bf21,
                                 file = "model/m_cohesion_an_bf21.rds")

6.2.1.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_cohesion_an <- dh_check_brms(m_cohesion_an, quantreg = TRUE)

## 21年より前
dh_cohesion_an_bf21 <- dh_check_brms(m_cohesion_an_bf21, quantreg = TRUE)


多重共線性の問題もありません。

## 全期間   
check_collinearity(m_cohesion_an)
## 21年より前
check_collinearity(m_cohesion_an_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_cohesion_an, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_cohesion_an_bf21, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_cohesion_an)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_cohesion_an_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


6.2.1.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。


続いて、2021年のデータを除いたモデルについて確認する。


6.2.1.5 結果の図示

最後に、結果の図示をします。

まず、TYの有無についてです。

emm_cohesion_an_TY <- emmeans(m_cohesion_an,
                              ~TY_presence)

cohesion_data_an %>% 
  ggplot(aes(x = TY_presence, y = log(mean_RG))) +
  geom_violin(bw = 0.2,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_cohesion_an_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_cohesion_an_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  geom_signif(comparisons = list(c("0", "1")),
              y_position = 2.5,
              annotations = "***") +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
  labs(x = "", y = "log(mean resting group size)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_cohesion_an_TY

p_cohesion_an_TY

ggsave("figures/p_cohesion_an_TY.tif", p_cohesion_an_TY,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


続いて、ITについてです。

emm_cohesion_an_IT <- emmeans(m_cohesion_an_bf21,
                              ~IT_presence)

cohesion_data_an_bf21 %>% 
  ggplot(aes(x = IT_presence, y = log(mean_RG))) +
  geom_violin(bw = 0.2,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_cohesion_an_IT %>% 
               data.frame() %>% 
               mutate(IT_presence = as.factor(IT_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_cohesion_an_IT %>% 
               data.frame() %>% 
               mutate(IT_presence = as.factor(IT_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*IT* absent", "*IT* present")) +
  labs(x = "", y = "log(mean resting group size)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_cohesion_an_IT

p_cohesion_an_IT

ggsave("figures/p_cohesion_an_IT.tif", p_cohesion_an_IT,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


6.2.2 非交尾期

6.2.2.1 データの加工

続いて、非交尾期について分析を行うためのデータを作成します。連続変数は標準化します。

  • 休息ポイント数が10以上の個体追跡セッションのみを使用
  • ハドリングのポイントは除去
  • 地上休息のみを含む
cohesion_data_nm <- focal_combined_prox_final %>% 
  mutate(rate_female = no_female/max_female) %>%
  filter(str_detect(study_period, "nm")) %>% 
  replace_na(list(hud = 0)) %>% 
  filter(hud != "1") %>% 
  filter(RG == 1 & T_G == "G") %>% 
  group_by(date, no_focal, subject, study_period, 
           no_female, TY_presence,
           no_kin, IT_presence) %>% 
  summarise(sum_RG = sum(RG_female_plus1),
            n = n()) %>% 
  ungroup() %>% 
  mutate(TY_presence = as.factor(TY_presence),
         IT_presence = as.factor(IT_presence)) %>%
  filter(study_period != "m18") %>% 
  filter(n >= 10) %>% 
  mutate(kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         mean_RG = sum_RG/n)

cohesion_data_nm


共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。

cohesion_data_nm %>% 
  pivot_longer(cols = c(no_female,  no_kin),
               names_to = "type",
               values_to = "value") %>% 
  ggplot(aes(x = value, y = mean_RG)) +
  geom_point(shape = 1) +
  facet_rep_wrap(~ type,
                 scales = "free") +
  theme_bw(base_size = 12) +
  theme(aspect.ratio = 1) +
  labs(x = "", y = "Mean resting group size")


6.2.2.2 モデリング

それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。

  • 応答変数: log(平均休息集団サイズ(メスのみ))
  • 分布: Studentのt分布
  • 説明変数: TYの在不在(TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
  • ランダム切片: subject
## 2019~2021年
m_cohesion_nm <- brm(log(mean_RG) ~ TY_presence + kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd"),
                                           prior(student_t(4,0,10), class = "nu"),
                                           prior(student_t(4,0,10), class = "sigma")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = cohesion_data_nm,
                                 file = "model/m_cohesion_nm.rds")

## 切断モデル  
m_cohesion_nm_trunc <- brm(log(mean_RG)|trunc(lb = 0) ~ TY_presence + 
                       kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 7000, warmup = 2000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd"),
                                           prior(student_t(4,0,10), class = "nu"),
                                           prior(student_t(4,0,10), class = "sigma")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = cohesion_data_nm,
                                 file = "model/m_cohesion_nm.rds")

## 2019~2020年
cohesion_data_nm_bf21 <- cohesion_data_nm %>% 
  filter(study_period != "m22") %>% 
  mutate(kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         mean_RG = sum_RG/n)

m_cohesion_nm_bf21 <- brm(log(mean_RG) ~ TY_presence + kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = cohesion_data_nm_bf21,
                                 file = "model/m_cohesion_nm_bf21.rds")

## 切断モデル  
m_cohesion_nm_bf21_trunc <- brm(log(mean_RG)|trunc(lb = 0) ~ TY_presence + IT_presence + ntm_std + est_std + 
                                  female_std + kin_std  + study_period + (1|subject),
                                family = student,
                                iter = 11000, warmup = 1000, seed = 123,
                                prior = c(prior(student_t(4,0,5), class = "b"),
                                          prior(student_t(4,0,10), class = "Intercept"),
                                          prior(student_t(4,0,10), class = "sd")),
                                control = list(adapt_delta = 0.99, max_treedepth = 20),
                                backend = "cmdstanr",
                                data = cohesion_data_nm_bf21,
                                file = "model/m_cohesion_nm_bf21.rds")

6.2.2.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_cohesion_nm <- dh_check_brms(m_cohesion_nm, quantreg = TRUE)

## 21年より前
dh_cohesion_nm_bf21 <- dh_check_brms(m_cohesion_nm_bf21, quantreg = TRUE)


多重共線性の問題もありません。

## 全期間   
check_collinearity(m_cohesion_nm)
## 21年より前
check_collinearity(m_cohesion_nm_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_cohesion_nm, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_cohesion_nm_bf21, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_cohesion_nm)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_cohesion_nm_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


6.2.2.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。


続いて、2021年のデータを除いたモデルについて確認する。


6.2.2.5 結果の図示

最後に、結果の図示をします。

まず、TYの有無についてです。

emm_cohesion_nm_TY <- emmeans(m_cohesion_nm,
                              ~TY_presence)

cohesion_data_nm %>% 
  ggplot(aes(x = TY_presence, y = log(mean_RG))) +
  geom_violin(bw = 0.2,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_cohesion_nm_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_cohesion_nm_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
  labs(x = "", y = "log(mean resting group size)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_cohesion_nm_TY

p_cohesion_nm_TY

ggsave("figures/p_cohesion_nm_TY.tif", p_cohesion_nm_TY,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


6.3 近接個体数の分析(休息中)

ここでは、TYやITの有無を含む様々な要因が個体追跡中の周辺メス数に影響するかを検討します。

6.3.1 交尾期(非発情メス)

6.3.1.1 データの加工

まず、非発情メスについて分析を行うためのデータを作成します。連続変数は標準化します。

  • 休息ポイント数が10以上の個体追跡セッションのみを使用
  • ハドリングのポイントは除去
  • 地上休息のみを含む
prox_data_an_RG <- focal_combined_prox_final %>% 
  mutate(rate_female = no_female/max_female) %>%
  filter(rs2 == 0) %>% 
  replace_na(list(hud = 0)) %>% 
  filter(hud != "1") %>% 
  filter(RG == 1 & T_G == "G") %>% 
  group_by(date, no_focal, subject, study_period, rate_agg_6yo, no_ntm, 
           no_est, no_male, no_female, TY_presence,
           no_kin, IT_presence, CSI_TY) %>% 
  summarise(sum_prox = sum(x3m_female + 1),
            n = n()) %>% 
  ungroup() %>% 
  mutate(TY_presence = as.factor(TY_presence),
         IT_presence = as.factor(IT_presence)) %>%
  filter(study_period != "m18") %>% 
  filter(n >= 10) %>% 
  mutate(ntm_std = standardize(no_ntm),
         male_std = standardize(no_male),
         est_std = standardize(no_est),
         kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         rate_agg_std = standardize(rate_agg_6yo),
         CSI_TY_std = standardize(CSI_TY),
         mean_prox = sum_prox/n)

prox_data_an_RG


共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。

prox_data_an_RG %>% 
  pivot_longer(cols = c(no_ntm, no_est, no_female, rate_agg_6yo, kin_std, CSI_TY),
               names_to = "type",
               values_to = "value") %>% 
  ggplot(aes(x = value, y = mean_prox)) +
  geom_point(shape = 1) +
  facet_rep_wrap(~ type,
                 scales = "free") +
  theme_bw(base_size = 12) +
  theme(aspect.ratio = 1) +
  labs(x = "", y = "Mean resting group size")


6.3.1.2 モデリング

それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。

  • 応答変数: log(平均休息集団サイズ(メスのみ))
  • 分布: Studentのt分布
  • 説明変数: TYの在不在(TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
  • ランダム切片: subject
## 2019~2021年
m_prox_an_RG <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std + 
                      female_std + kin_std  + study_period +  (1|subject) + (1|date),
                    family = student,
                    iter = 11000, warmup = 1000, seed = 123,
                    prior = c(prior(student_t(4,0,5), class = "b"),
                              prior(student_t(4,0,10), class = "Intercept"),
                              prior(student_t(4,0,10), class = "sd"),
                              prior(student_t(4,0,10), class = "nu"),
                              prior(student_t(4,0,10), class = "sigma")),
                    control = list(adapt_delta = 0.99, max_treedepth = 20),
                    backend = "cmdstanr",
                    data = prox_data_an_RG,
                    file = "model/m_prox_an_RG.rds")

## 2019~2020年
prox_data_an_RG_bf21 <- prox_data_an_RG %>% 
  filter(study_period != "m21") %>% 
  mutate(ntm_std = standardize(no_ntm),
         male_std = standardize(no_male),
         est_std = standardize(no_est),
         kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         rate_agg_std = standardize(rate_agg_6yo),
         CSI_TY_std = standardize(CSI_TY),
         mean_prox = sum_prox/n)

m_prox_an_RG_bf21 <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std + 
                       female_std + kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = prox_data_an_RG_bf21,
                                 file = "model/m_prox_an_RG_bf21.rds")

6.3.1.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_prox_an_RG <- dh_check_brms(m_prox_an_RG, quantreg = TRUE)

## 21年より前
dh_prox_an_RG_bf21 <- dh_check_brms(m_prox_an_RG_bf21, quantreg = TRUE)


多重共線性の問題もありません。

## 全期間   
check_collinearity(m_prox_an_RG)
## 21年より前
check_collinearity(m_prox_an_RG_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_prox_an_RG, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_prox_an_RG_bf21, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_prox_an_RG)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_prox_an_RG_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


6.3.1.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。


続いて、2021年のデータを除いたモデルについて確認する。


6.3.1.5 結果の図示

最後に、結果の図示をします。

まず、TYの有無についてです。

emm_prox_an_RG_TY <- emmeans(m_prox_an_RG,
                          ~TY_presence)

prox_data_an_RG %>% 
  ggplot(aes(x = TY_presence, y = log(mean_prox))) +
  geom_violin(bw = 0.2,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_prox_an_RG_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_prox_an_RG_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  geom_signif(comparisons = list(c("0", "1")),
              y_position = 2.1,
              annotations = "***") +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
  labs(x = "", y = "log(mean no. of females within 3m)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_prox_an_RG_TY

p_prox_an_RG_TY

ggsave("figures/p_prox_an_RG_TY.tif", p_prox_an_RG_TY,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


続いて、ITについてです。

emm_prox_an_RG_IT <- emmeans(m_prox_an_RG_bf21,
                             ~IT_presence)

prox_data_an_RG_bf21 %>% 
  ggplot(aes(x = IT_presence, y = log(mean_prox))) +
  geom_violin(bw = 0.2,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_prox_an_RG_IT %>% 
               data.frame() %>% 
               mutate(IT_presence = as.factor(IT_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_prox_an_RG_IT %>% 
               data.frame() %>% 
               mutate(IT_presence = as.factor(IT_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*IT* absent", "*IT* present")) +
  labs(x = "", y = "log(mean resting group size)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_prox_an_RG_IT

p_prox_an_RG_IT

ggsave("figures/p_prox_an_RG_IT.tif", p_prox_an_RG_IT,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


6.3.2 非交尾期

6.3.2.1 データの加工

続いて、非交尾期について分析を行うためのデータを作成します。連続変数は標準化します。

  • 休息ポイント数が10以上の個体追跡セッションのみを使用
  • ハドリングのポイントは除去
  • 地上休息のみを含む
prox_data_nm_RG <- focal_combined_prox_final %>% 
  mutate(rate_female = no_female/max_female) %>%
  filter(str_detect(study_period, "nm")) %>% 
  replace_na(list(hud = 0)) %>% 
  filter(hud != "1") %>% 
  filter(RG == 1 & T_G == "G") %>% 
  group_by(date, no_focal, subject, study_period, 
           no_female, TY_presence,
           no_kin, IT_presence) %>% 
  summarise(sum_prox = sum(x3m_female + 1),
            n = n()) %>% 
  ungroup() %>% 
  mutate(TY_presence = as.factor(TY_presence),
         IT_presence = as.factor(IT_presence)) %>%
  filter(study_period != "m18") %>% 
  filter(n >= 10) %>% 
  mutate(kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         mean_prox = sum_prox/n)

prox_data_nm_RG


共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。

prox_data_nm_RG %>% 
  pivot_longer(cols = c(no_female,  no_kin),
               names_to = "type",
               values_to = "value") %>% 
  ggplot(aes(x = value, y = mean_prox)) +
  geom_point(shape = 1) +
  facet_rep_wrap(~ type,
                 scales = "free") +
  theme_bw(base_size = 12) +
  theme(aspect.ratio = 1) +
  labs(x = "", y = "Mean resting group size")


6.3.2.2 モデリング

それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。

  • 応答変数: log(平均休息集団サイズ(メスのみ))
  • 分布: Studentのt分布
  • 説明変数: TYの在不在(TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
  • ランダム切片: subject
## 2019~2021年
m_prox_nm_RG <- brm(log(mean_prox) ~ TY_presence + kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd"),
                                           prior(student_t(4,0,10), class = "nu"),
                                           prior(student_t(4,0,10), class = "sigma")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = prox_data_nm_RG,
                                 file = "model/m_prox_nm_RG.rds")

6.3.2.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_prox_nm_RG <- dh_check_brms(m_prox_nm_RG, quantreg = TRUE)


多重共線性の問題もありません。

## 全期間   
check_collinearity(m_prox_nm_RG)

bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_prox_nm_RG, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_prox_nm_RG)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


6.3.2.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。


6.3.2.5 結果の図示

最後に、結果の図示をします。

まず、TYの有無についてです。

emm_prox_nm_RG_TY <- emmeans(m_prox_nm_RG,
                          ~TY_presence)

prox_data_nm_RG %>% 
  ggplot(aes(x = TY_presence, y = log(mean_prox))) +
  geom_violin(bw = 0.2,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_prox_nm_RG_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_prox_nm_RG_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
  labs(x = "", y = "log(mean resting group size)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_prox_nm_RG_TY

p_prox_nm_RG_TY

ggsave("figures/p_prox_nm_RG_TY.tif", p_prox_nm_RG_TY,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


6.4 近接個体数の分析(採食中)

最後に、TYやITの有無を含む様々な要因が個体追跡中の周辺メス数に影響するかを検討します。

6.4.1 交尾期(非発情メス)

6.4.1.1 データの加工

まず、非発情メスについて分析を行うためのデータを作成します。連続変数は標準化します。

  • 休息ポイント数が10以上の個体追跡セッションのみを使用
  • ハドリングのポイントは除去
  • 地上休息のみを含む
prox_data_an_F <- focal_combined_prox_final %>% 
  mutate(rate_female = no_female/max_female) %>%
  filter(rs2 == 0) %>% 
  replace_na(list(hud = 0)) %>% 
  filter(hud != "1") %>% 
  filter(activity == "F" & T_G == "G") %>% 
  group_by(date, no_focal, subject, study_period, rate_agg_6yo, no_ntm, 
           no_est, no_male, no_female, TY_presence,
           no_kin, IT_presence, CSI_TY) %>% 
  summarise(sum_prox = sum(x3m_female + 1),
            n = n()) %>% 
  ungroup() %>% 
  mutate(TY_presence = as.factor(TY_presence),
         IT_presence = as.factor(IT_presence)) %>%
  filter(study_period != "m18") %>% 
  filter(n >= 10) %>% 
  mutate(ntm_std = standardize(no_ntm),
         male_std = standardize(no_male),
         est_std = standardize(no_est),
         kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         rate_agg_std = standardize(rate_agg_6yo),
         CSI_TY_std = standardize(CSI_TY),
         mean_prox = sum_prox/n)

prox_data_an_F


共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。

prox_data_an_F %>% 
  pivot_longer(cols = c(no_ntm, no_est, no_female, rate_agg_6yo, kin_std, CSI_TY),
               names_to = "type",
               values_to = "value") %>% 
  ggplot(aes(x = value, y = mean_prox)) +
  geom_point(shape = 1) +
  facet_rep_wrap(~ type,
                 scales = "free") +
  theme_bw(base_size = 12) +
  theme(aspect.ratio = 1) +
  labs(x = "", y = "Mean resting group size")


6.4.1.2 モデリング

それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。

  • 応答変数: log(平均休息集団サイズ(メスのみ))
  • 分布: Studentのt分布
  • 説明変数: TYの在不在(TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
  • ランダム切片: subject
## 2019~2021年
m_prox_an_F <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std + 
                      female_std + kin_std  + study_period +  (1|subject) + (1|date),
                    family = student,
                    iter = 11000, warmup = 1000, seed = 123,
                    prior = c(prior(student_t(4,0,5), class = "b"),
                              prior(student_t(4,0,10), class = "Intercept"),
                              prior(student_t(4,0,10), class = "sd"),
                              prior(student_t(4,0,10), class = "nu"),
                              prior(student_t(4,0,10), class = "sigma")),
                    control = list(adapt_delta = 0.99, max_treedepth = 20),
                    backend = "cmdstanr",
                    data = prox_data_an_F,
                    file = "model/m_prox_an_F.rds")

## 2019~2020年
prox_data_an_F_bf21 <- prox_data_an_F %>% 
  filter(study_period != "m21") %>% 
  mutate(ntm_std = standardize(no_ntm),
         male_std = standardize(no_male),
         est_std = standardize(no_est),
         kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         rate_agg_std = standardize(rate_agg_6yo),
         CSI_TY_std = standardize(CSI_TY),
         mean_prox = sum_prox/n)

m_prox_an_F_bf21 <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std + 
                       female_std + kin_std  + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 123,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = prox_data_an_F_bf21,
                                 file = "model/m_prox_an_F_bf21.rds")

6.4.1.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_prox_an_F <- dh_check_brms(m_prox_an_F, quantreg = TRUE)

## 21年より前
dh_prox_an_F_bf21 <- dh_check_brms(m_prox_an_F_bf21, quantreg = TRUE)


多重共線性の問題もありません。

## 全期間   
check_collinearity(m_prox_an_F)
## 21年より前
check_collinearity(m_prox_an_F_bf21)

bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_prox_an_F, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)

## 21年より前
pp_check(m_prox_an_F_bf21, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_prox_an_F)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)

data.frame(Rhat = brms::rhat(m_prox_an_F_bf21)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


6.4.1.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。


続いて、2021年のデータを除いたモデルについて確認する。


6.4.1.5 結果の図示

最後に、結果の図示をします。

まず、TYの有無についてです。

emm_prox_an_F_TY <- emmeans(m_prox_an_F,
                          ~TY_presence)

prox_data_an_F %>% 
  ggplot(aes(x = TY_presence, y = log(mean_prox))) +
  geom_violin(bw = 0.1,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_prox_an_F_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_prox_an_F_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  geom_signif(comparisons = list(c("0", "1")),
              y_position = 1.3,
              annotations = "***") +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
  labs(x = "", y = "log(mean no. of females within 3m)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_prox_an_F_TY

p_prox_an_F_TY

ggsave("figures/p_prox_an_F_TY.tif", p_prox_an_F_TY,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


続いて、ITについてです。

emm_prox_an_F_IT <- emmeans(m_prox_an_F_bf21,
                             ~IT_presence)

prox_data_an_F_bf21 %>% 
  ggplot(aes(x = IT_presence, y = log(mean_prox))) +
  geom_violin(bw = 0.07,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_prox_an_F_IT %>% 
               data.frame() %>% 
               mutate(IT_presence = as.factor(IT_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_prox_an_F_IT %>% 
               data.frame() %>% 
               mutate(IT_presence = as.factor(IT_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*IT* absent", "*IT* present")) +
  labs(x = "", y = "log(mean resting group size)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_prox_an_F_IT

p_prox_an_F_IT

ggsave("figures/p_prox_an_F_IT.tif", p_prox_an_F_IT,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


6.4.2 非交尾期

6.4.2.1 データの加工

続いて、非交尾期について分析を行うためのデータを作成します。連続変数は標準化します。

  • 休息ポイント数が10以上の個体追跡セッションのみを使用
  • ハドリングのポイントは除去
  • 地上休息のみを含む
prox_data_nm_F <- focal_combined_prox_final %>% 
  mutate(rate_female = no_female/max_female) %>%
  filter(str_detect(study_period, "nm")) %>% 
  replace_na(list(hud = 0)) %>% 
  filter(hud != "1") %>% 
  filter(activity == "F" & T_G == "G") %>% 
  group_by(date, no_focal, subject, study_period, 
           no_female, TY_presence,
           no_kin, IT_presence) %>% 
  summarise(sum_prox = sum(x3m_female + 1),
            n = n()) %>% 
  ungroup() %>% 
  mutate(TY_presence = as.factor(TY_presence),
         IT_presence = as.factor(IT_presence)) %>%
  filter(study_period != "m18") %>% 
  filter(n >= 10) %>% 
  mutate(kin_std = standardize(no_kin),
         female_std = standardize(no_female),
         mean_prox = sum_prox/n)

prox_data_nm_F


共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。

prox_data_nm_F %>% 
  pivot_longer(cols = c(no_female,  no_kin),
               names_to = "type",
               values_to = "value") %>% 
  ggplot(aes(x = value, y = mean_prox)) +
  geom_point(shape = 1) +
  facet_rep_wrap(~ type,
                 scales = "free") +
  theme_bw(base_size = 12) +
  theme(aspect.ratio = 1) +
  labs(x = "", y = "Mean resting group size")


6.4.2.2 モデリング

それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。

  • 応答変数: log(平均休息集団サイズ(メスのみ))
  • 分布: Studentのt分布
  • 説明変数: TYの在不在(TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
  • ランダム切片: subject
## 2019~2021年
m_prox_nm_F <- brm(log(mean_prox) ~ TY_presence + kin_std + study_period + (1|subject),
                                 family = student,
                                 iter = 11000, warmup = 1000, seed = 195,
                                 prior = c(prior(student_t(4,0,5), class = "b"),
                                           prior(student_t(4,0,10), class = "Intercept"),
                                           prior(student_t(4,0,10), class = "sd"),
                                           prior(student_t(4,0,10), class = "nu"),
                                           prior(student_t(4,0,10), class = "sigma")),
                                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                                 backend = "cmdstanr",
                                 data = prox_data_nm_F,
                                 file = "model/m_prox_nm_F.rds")

6.4.2.3 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。

## 全期間    
dh_prox_nm_F <- dh_check_brms(m_prox_nm_F, quantreg = TRUE)


多重共線性の問題もありません。

## 全期間   
check_collinearity(m_prox_nm_F)

bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_prox_nm_F, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_prox_nm_F)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


6.4.2.4 結果の確認

まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。


6.4.2.5 結果の図示

最後に、結果の図示をします。

まず、TYの有無についてです。

emm_prox_nm_F_TY <- emmeans(m_prox_nm_F,
                          ~TY_presence)

prox_data_nm_F %>% 
  ggplot(aes(x = TY_presence, y = log(mean_prox))) +
  geom_violin(bw = 0.2,
              scale = "width") +
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_prox_nm_F_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_prox_nm_F_TY %>% 
               data.frame() %>% 
               mutate(TY_presence = as.factor(TY_presence)),
             aes(y = emmean, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
  labs(x = "", y = "log(mean resting group size)") +
  theme(text = element_text(family = "Times New Roman"),
        axis.text.x = element_markdown(),
        aspect.ratio = 1) -> p_prox_nm_F_TY

p_prox_nm_F_TY

ggsave("figures/p_prox_nm_F_TY.tif", p_prox_nm_F_TY,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


References

Gabry, J., & Mahr, T. (2022). Bayesplot: Plotting for bayesian models. https://mc-stan.org/bayesplot
Hartig, F. (2022). DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. https://CRAN.R-project.org/package=DHARMa
Rodríguez-Sánchez, F. (2023). DHARMa.helpers: Helper functions to check models not (yet) directly supported by DHARMa.